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We report on simulations of reversible random sequen- 
tial adsorption of dimers on three different lattices: a one- 
dimensional lattice, a two-dimensional triangular lattice, and 
a two-dimensional triangular lattice with the nearest neigh- 
bors excluded. In addition to the adsorption of particles at 
a rate K + , we allow particles to leave the surface at a rate 
K~ . The results from the one-dimensional lattice model agree 
with previous results for the continuous parking lot model. 
In particular, the long-time behavior is dominated by collec- 
tive events involving two particles. We were able to directly 
confirm the importance of two-particle events in the simple 
two-dimensional triangular lattice. For the two-dimensional 
triangular lattice with the nearest neighbors excluded, the 
observed dynamics are consistent with this picture. The two- 
dimensional simulations were motivated by measurements of 
Ca ++ binding to Langmuir monolayers. The two cases were 
chosen to model the effects of changing pH in the experimen- 
tal system. 
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I. INTRODUCTION 

A large number of nonequilibrium systems can be qual- 
itatively described as a flux of particles impinging on a 
surface or line. Two heavily studied models of such sys- 
tems treat the particles as either fixed in place upon im- 
pact (random sequential adsorption) or as free to diffuse 
along the surface or line (random cooperative adsorption) 
|jj . One can also consider the deposition of particles that 
are free to desorb |^-^| ■ Some examples of the wide range 
of applicability of these models include coating problems, 
chemisorption, physisorption, the reaction of molecular 
species on surfaces and at interfaces, and the binding of 
ligands on polymer chains. Jamming is one of the com- 
mon occurrences in these systems that random sequential 
adsorption models effectively describe. Loosely speak- 
ing, a jammed system is one that is locked into a state 
of partial coverage because of adsorbate size or shape. 
In addition to the various adsorption processes, jamming 
occurs in a wide range of nonequilibrium situations, in- 
cluding glasses, granular materials, and traffic flow |i] ||. 
In spite of significant progress, no general framework ex- 
ists for the description of jamming phenomena. 

A particular realization of random sequential adsorp- 
tion is the parking lot model Jl|j9|-|l^] . In the irreversible 



version of this model, identical particles (cars) adsorb 
on a line (curb) at a rate K + . In this model, the phe- 
nomenon of jamming has been known for some time ||. 
A certain number of the parked cars leave a space that is 
too small to fit another car. These are referred to as bad 
parkers. The result is a density of cars along the curb 
that is less than one. The density of cars reached in the 
irreversible model is the jamming limit. 

In the reversible version, identical particles (cars) ad- 
sorb on a line (curb) at a rate K + and leave the line 
(curb) at a rate K~ . The removal of cars allows for 
adjustments in the bad parkers that relieve the jam- 
ming. Recently, there has been renewed interest in the re- 
versible case because of its successful application to com- 
paction in granular materials when generalized to three 
dimensions In this version, the "parking spots" are 
voids in the material that can be filled with particles. The 
dynamics of the reversible parking lot model for large val- 
ues of K = K + /K~ has a number of interesting features. 
Perhaps the most dramatic feature is the existence of two 
very different time scales for the evolution of the coverage 
fraction of particles fll4|| . First, there is a rapid approach 
to a coverage fraction that is equal to the jamming limit. 
This is followed by a slow relaxation to a larger steady 
state value. The slow relaxation is understood in terms of 
collective parking/leaving events involving multiple cars 

In this paper, we present the results for simulations 
of the reversible adsorption of dimers on: (1) a one- 
dimensional lattice, (2) a two-dimensional triangular lat- 
tice, and (3) a two-dimensional triangular lattice with the 
nearest neighbors excluded. The one-dimensional lattice 
model |p^ , p5| -^7| was chosen as a test case, and the re- 
sults are in good agreement with existing data. In par- 
ticular, our simulations confirm the importance of col- 
lective parking events in controlling the slow dynamics, 
as seen in Ref. [jlij. The two triangular lattice models 
exhibit differences in their time evolution that can be at- 
tributed to effects of bond orientation and packing on 
the collective events. The case without nearest-neighbor 
exclusion corresponds to attempting to cover the plane 
with a shape formed by two regular hexagons sharing a 
side. The nearest-neighbor excluded case corresponds to 
a tiling of distorted hexagons that cover multiple sites. 

The rest of the paper is organized as follows. Section 
II describes the details of the simulations. Section III 
presents the results for the one-dimensional model. Sec- 
tion IV presents the results for the two triangular lattices. 
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The simulations were motivated in part by experimental 
measurements of the viscosity of Langmuir monolayers. 
A brief description of the experimental system and its 
relationship to the simulations presented here is given in 
Section V. The results are discussed and summarized in 
Section VI. 



II. SIMULATION DETAILS 

For the one-dimensional simulations, a line of 32000 
particles was used. Both of the triangular lattices con- 
sisted of a grid of 1000 x 1000 particles. To distin- 
guish between the two-dimensional models, we intro- 
duce the following nomenclature. Model A will refer to 
the triangular lattice without nearest-neighbor exclusion. 
Model B will refer to the triangular lattice with nearest- 
neighbors excluded from binding. Particles are taken to 
bind to two neighboring sites on the lattice, forming a 
dimer. The binding occurs at a rate K + , and particles 
leave the surface at a rate of K~ . 

At each step in the simulation, a site was chosen at 
random. Then, a random number between and 1 was 
compared with the ratio K + /(K + + K~) to determine 
whether a binding or unbinding event was attempted. 
For unbinding events, if the chosen site was part of a 
bond, the bond was broken; otherwise, no action took 
place. For binding events, a nearest neighbor was ran- 
domly selected. A binding event occurred only if both 
sites were allowed binding sites. The definition of al- 
lowed binding site depends on the model. For the one- 
dimensional and Model A cases, an allowed site is any 
site that is not part of a bond. For Model B, if either site 
is part of a bond or the nearest neighbor of a bound site, 
binding is not allowed. It is important to note that the 
number of new bonds created is directly proportional to 
the number of allowed sites, which is not the same as the 
number of open sites. The number of desorption events 
is still directly proportional to the coverage fraction. The 
coverage fraction, p, is defined as the ratio of sites that 
are part of a bond to the total number of sites. 

A schematic of each of the model systems with ex- 
amples of bound sites is shown in Fig. 1. It is impor- 
tant to notice the different spatial structures in Model 
A and Model B. In Model A, complete coverage corre- 
sponds to all sites being part of a bond. In Model B, 
perfect coverage of the system corresponds to a tile of 
distorted hexagons that are composed of both empty and 
bound sites. This results in a maximal coverage frac- 
tion p m ax = 0.4. For both the one-dimensional case and 
Model A, 

Pmax — 1-0. In this paper, p(oo) will desig- 
nate the steady state value of the fractional coverage, 
and pjam will refer to p(oo) in the case K~ = 0, i.e. the 
jamming limit. 



(a) 



(b) 




(c) 




FIG. 1. (a) Example of bonds in the one-dimensional lat- 
tice model. The sites are represented by the dots, and the 
incoming particles are represented by solid lines. The parti- 
cles form a bond between two neighboring sites, (b) Example 
of bonds in Model A. In this case, particles can bind any two 
nearest-neighbor sites that are not already part of a bond, (c) 
Example of bonds in Model B. In this case, nearest-neighbors 
of a bound site are not allowed to form bonds. Examples of 
such sites are represented by open circles. 



III. ONE-DIMENSIONAL SIMULATION 

Figure 2 shows p as a function of iteration step for 
selected values of K = K + /K~. We include the case 
K~ = 0, which gives Pj am = 0.86474. For comparison, 
analytic calculations give Pj am = 0.86466 p^jl Tj ], The 
original work on the reversible parking lot modelp2| pro- 
posed a mean-field description of the dynamics that can 
be expressed in terms of the average density, or fractional 
surface coverage, p. Both the continuous and lattice ver- 
sions of the parking lot model were considered. Figure 
3 shows the steady state value of p for the values of K 
plotted in Fig. 2. The solid curve in Fig. 3 is the value 
for p(oo) for dimers binding to a one-dimensional lattice, 
as determined by the following equation from Ref. |12]] : 



p{oo) = 1 - (K-/K + y/ 2 /2. 



(1) 
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The agreement between our simulations and Eq. 1 con- 
firms the mean field prediction for the equilibrium values 
of p. However, as with the continuous parking lot model 
||l4|| , the mean-field description is unable to accurately 
predict the time evolution of p. This can be seen in Fig. 
2 where the two time scales controlling the evolution of 
p are evident for K > 10. The system rapidly reaches 
Pjam, and then slowly approaches its equilibrium value. 
As K goes to infinity, p(oo) approaches one, but the time 
to reach equilibrium approaches infinity. This is in agree- 
ment with results for the continuous parking lot model 
reported in Ref. |l4[. 
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FIG. 2. Shown here is the coverage fraction p as a func- 
tion of the number of iterations for the one-dimensional 
model. The dashed line represents the evolution for K~ — 0. 
The other curves are for: solid line (K = 5), open trian- 
gles (K = 40), open squares (K — 200), and open circles 
(K = 1000). 

We have found that the explanation of the two distinct 
time scales reported in Ref. [Q applies to the discrete 
case as well. Essentially, collective events are responsi- 
ble for the evolution of p for p > pj am - 

In Ref. [14 , 

the authors calculated the transition rates for two good 
particles to one bad particle and one bad particle to two 
good particles and found that these rates account for the 
additional slow time scales. In contrast, we directly mon- 
itor the transitions as part of the simulation. The reason 
such transitions result in an additional slow time scale 
can be understood in terms of the following argument. 
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FIG. 3. The solid line is the mean-field prediction for p(oo) 
as a function of K. The symbols are the values of p(oo) taken 
from Fig. 2. 

As discussed in the introduction, when K~ = 0, jam- 
ming occurs because of "bad parkers" that leave empty 
space. For the one-dimensional lattice, empty space 
refers to a single site that is unable to bond. An ex- 
ample is shown in Fig. 4a. For small values of K~ , bad 
parkers initially occur at essentially the same rate as for 
K~ = because very few particles desorb. Therefore, 
the coverage fraction for the system quickly approaches 
a value of pj am - Even when a value of pj am is reached, 
the rare desorption event is generally followed immedi- 
ately by a readsorption because K + is so large. The 
total number of particles is not changed by these events. 
However, when one bad parker desorbs and two particles 
adsorb in the opened good locations, then the number 
of particles is increased by one. Likewise, if two good 
parkers unbind and one bad parker binds, the number 
of particles is decreased by one. Because these events in- 
volve multiple particle transitions, they occur on a longer 
time scale then simple adsorption/desorption events. 

For the one-dimensional discrete case, one can identify 
the relevant good to bad and bad to good transition that 
involve only two good parkers. These are illustrated in 
Figs. 4b and Fig. 4c. As these events are expected 
to dominate the dynamics, one can write the following 
equation for the evolution of p once the jamming limit 
has been reached: 



dp/dt — Rb g — Rgb + h.o.t. 



(2) 



Here Rb g and R g b are the rates of bad to good and good 
to bad transitions respectively, and h.o.t are collective 
transitions that a larger number of particles, and hence, 
occur at a slower rate. 
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FIG. 4. (a) This illustrates the concept of a "bad" parker. 
The bond at location A is the bad parker, as it leaves two sites 
free, (b) This illustrates a transition from one bad parker to 
two good parkers. The dashed line represents the original bad 
parker at location A that desorbs. Then, the two good park- 
ers, B and C, adsorb, (c) This illustrates the two good parkers 
to one bad parker transition. In this case, the two parkers at 
B and C leave and one attaches at A. (d) This illustrates a 
spatial arrangement that corresponds to a two-to-three parti- 
cle transition. The two possible bond distributions are given 
by the dashed lines at sites A, B, and C (the good parkers) 
and the solid lines at sites D and E, respectively. 

We were able to track the bad to good and good to bad 
transitions during the simulation. This was accomplished 
by converting the particle sites to an array of bond lo- 
cations. Each location between two sites was assigned a 
value of 1 if a bond was present and is there was no 
bond. For example, the solid lines in Fig. 4b would be 
represented by the string 1010101. Notice, by definition, 
between any two bonds there is an open space, so the 
completely filled system is represented by 1010101010.... 
The string of bond locations was saved at step i and i+ A. 
Each bond location was taken as the initial digit in a 
seven digit string, and these strings were compared for 
steps i and i + A. We counted the following transitions: 

1010101 1001001. 

These transitions correspond to two good to one bad and 
one bad to two good, as discussed in Fig. 4c and Fig. 4b 
respectively. The choice of A is important. If A is too 
small, the transitions do not have enough time to com- 
plete. For example, in the extreme limit of choosing A 
to be a single time step, it is not possible to have multi- 
particle events, but the total number of bound sites can 
change by one. Essentially, A must be large enough for 
the multiparticle transitions to have time to complete. 
For A large enough, the recorded number of transitions 
is essentially independent of A. For the data reported 



here, we used A = 2 x 10 6 . 

In addition to counting multiparticle transition, we 
also recorded the total change in p. Figure 5 compares 
the actual value of p as a function of the number of it- 
erations with the value obtained using Eq. 2 and the 
computed number of bad to good and good to bad tran- 
sitions. Once the jamming limit is reached, the bad to 
good and good to bad transitions account for 94.3% of 
the change in p, confirming the general idea behind Eq. 
2. An additional 3.2% of the change in p is accounted 
for by considering a single class of three particle tran- 
sitions where three good parkers were replaced by two 
parkers, and the reverse process. These were counted by 
considering 9 digit strings and looking for the transition: 

101010101 ^> 100101001. 

This curve is also plotted in Fig. 5. The spatial ar- 
rangement corresponding to this transition is shown in 
Fig. 4d. It is important to note that when bonds exist 
at site D and E, the only way to increase the number of 
bound sites in this region is for two particles to desorb 
and three particles adsorb at sites A, B and C. 
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FIG. 5. Shown here is the coverage fraction p as a function 
of the number of iterations for the one-dimensional model 
and K — 1000 (solid curve). The curve is plotted starting at 
the end of the jamming limit plateau. Also plotted are two 
curves that are obtained by numerically integrating Eq. 2 
using the coverage fraction at the jamming limit as the initial 
state. The dashed line is the result when only the rates for 
the good to bad and bad to good transitions are included 
in Eq. 2. These transitions are described in Fig. 4. The 
dotted line shows the improvement at late times by including 
a single higher-order transition involving three good parkers 
converting to two parkers. 
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IV. TWO-DIMENSIONAL SIMULATIONS 



The results of the simulation for p as a function of 
iteration step for the adsorption of dimers on a two- 
dimensional triangular lattice (Model A) and on a two- 
dimensional triangular lattice with nearest-neighbor ex- 
clusion (Model B) are presented in Figs. 6 and 7, respec- 
tively. For Model A, previous simulations have found a 
jamming limit of 0.9243 JL8| . Our simulations give a value 
of 0.9120. For Model B, we find a jamming limit of 0.275. 
Recall that complete coverage in this case corresponds to 
p = 0.4. We are not aware of any previous work on a 
Model B type simulation. However, by appropriately in- 
cluding the empty nearest-neighbor sites in the definition 
of p, we can compare to simulations involving n-mers of 
length 6 that cover a hexagonal patch. These simulations 
find a jamming limit of 0.6847, and our converted value 
is 0.6875 g§. 
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FIG. 6. Shown here is the coverage fraction p as a func- 
tion of the number of iterations for Model A. The solid line 
represents the evolution for K~ =0. The other curves are 
for: solid diamonds (K = 200), crosses [K = 500), open 
squares (K = 1000), open circles (K = 2000), open triangles 
(K = 5000), and plus signs (K = 10000). 

The results for the two-dimensional cases are qualita- 
tively similar to the one-dimensional case. One observes 
multiple time scales: a rapid approach to the jamming 
limit and a slow relaxation to the steady-state value. This 
suggests that the same picture of multiparticle transitions 
will apply to the two-dimensional system. However, in 
contrast to the one-dimensional case, the identification 
of collective transitions is significantly more complex for 
Model A and B because of the number of arrangements 
due to differing orientations of the bonds that can pro- 
duce bad parkers. However, we did carry out a limited 
analysis for the case of Model A. 
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FIG. 7. Shown here is the coverage fraction p as a func- 
tion of the number of iterations for Model B. The dashed 
line represents the evolution for K~ = 0. The other curves 
are for: crosses (K — 200), open triangles (K — 500), solid 
line (K — 1000), closed circles (K = 2000), closed squares 
(K = 5000), and plus signs (K = 10000). 

The method used to track multiparticle events in 
Model A was similar in concept to the one dimensional 
case. However, because the bonds have orientation, we 
compared the actual sites instead of the bonds for three 
classes of transitions: 

0110 1111 

1 11 
1 1 1 

10 11 
1 1 1 

In this case, occupied sites are represented by 1 and un- 
occupied sites are represented by 0. Using sites instead 
of bonds results in some differences between the methods 
used in the two-dimensional and one-dimensional cases. 
First, the transitions counted in this manner correspond 
to classes of transitions in the following sense. Because 
we track sites and not bonds, two neareset neighbor sites 
can be occupied either because they share a bond or be- 
cause of two neighboring bonds that are at an angle to 
the line being considered. So, the first class of transi- 
tions includes the transitions that are exactly analogous 
to the one-dimensional good to bad transitions. But, it 
also includes multiparticle transitions that involve bonds 
at an angle to the horizontal and that successfully fill 
the empty sites along the horizontal. Second, the offset 
of the l's and 0's in the second two classes of transi- 
tions are important and reflect the underlying hexagonal 
lattice. Note that because only nearest-neighbor bond- 
ing is allowed, the diagonal connecting the two zeros in 
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each case is not an allowed binding site. Finally, for the 
two-dimensional case, we exploited the hexagonal sym- 
metry of the problem, and multiplied the rate for the 
first type of transition by three. The rates of the sec- 
ond two transitions are multiplied by 3/2 to account for 
double counting. The value of A was chosen in a similar 
fashion to the one-dimensional case and corresponded to 
a constant interval in Ap = 0.005. 

The results for p as a function of the number of itera- 
tions and the value of p computed from Eq. 2 using just 
these three classes of events defined above are plotted in 
Fig. 8. One striking feature of Fig. 8 is the fact that 
the two-particle events we identified account for nearly 
100% of the dynamics until the number of steps reaches 
approximately 2 x 10 9 . At this point, the coverage con- 
tinues to grow, only there is essentially no change due to 
the identified two-particle events. This strongly suggests 
other two-particle events or higher-order events involving 
more than two particles are becoming important. 
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FIG. 8. In this figure, the curves are the coverage fraction p 
versus number of iterations for Model A. A reduced grid of 600 
by 600 was used to facilitate the counting of collective events. 
The dashed line is the result of the simulation for K — 5000. 
The solid line is the result obtained by tracking two-particle 
events and integrating Eq. 2 starting at the jamming limit. 
The symbols are plotted against the right-hand axis and give 
the ratio of the change in p as computed by the two methods. 
The agreement between the two methods is excellent until 
approximately 2 x 10 9 steps. At this point, the contribution 
to the dynamics of two-particle events decreases dramatically. 



V. POSSIBLE APPLICATION TO LANGMUIR 
MONOLAYERS 

An obvious question is: do the triangular lattices con- 
sidered here apply to any experimental systems? There is 
indirect evidence that the models discussed here are rel- 



evant to the binding of Ca ++ ions to a Langmuir mono- 
layer. Langmuir monolayers are composed of insoluble, 
amphiphilic molecules that are confined to the air-water 
interface JT9). They exhibit the usual gas, liquid, and 
solid phases, as well as a large number of two-dimensional 
analogs of smectic phases ]2(J . Many of these phases are 
hexatic, with the molecules locally arranged on a dis- 
torted hexagonal lattice. When Ca ++ is present in the 
water, it can bind two fatty-acid molecules together. This 
substantially alters a number of the physical properties of 
the monolayer, such as the lattice spacing and the viscos- 
ity H]. Existing measurements (2^] and models (23| of 
Ca + ~binding have focused on the equilibrium coverage 
fraction. However, the measurements have focused on 
time scales of one hour or less. The coverage fraction de- 
pends strongly on pH, which is understandable in terms 
of the degree of ionization of the fatty acid headgroup. 
At low values of the pH, essentially all of the fatty acid 
molecules are neutral, and the Ca ++ ions do not bind. 
As the pH is increased, an increasing number of fatty 
acid molecules become charged, and the Ca ++ ions are 
free to bind to the monolayer. 

The possible relevance of Model A and B to the fatty 
acid monolayers is based on viscosity measurements as a 
function of time in the presence of Ca ++ for the hexatic 
phase of a particular fatty acid J24L Figure 9 reproduces 
one set of data from Fig. 2 of Ref. |24| , illustrating a typi- 
cal time evolution of the viscosity. The viscosity increases 
2 orders of magnitude over 15 hours. The time evolution 
can be divided into three distinct regions: an initial rapid 
rise in viscosity within the first hour, a slower rise in vis- 
cosity covering 5-6 hours, and a final even slower rise 
in viscosity. For comparison, the computed fractional 
coverage of p is shown in Fig. 9 versus the number of it- 
erations. In this case, we have used a linear scale for the 
number of iterations. The previous plots all used a log- 
arithmic scale. The time evolution of the Ca ++ binding 
exhibits the three general regions present in the viscosity 
data, and as such, provides a natural explanation for the 
effect. 

There are a number of points with regard to the con- 
nections between the model and the monolayer experi- 
ments. The simulations are consistent with the fact that 
previous measurements of Ca ++ binding do not observe 
multiple time scales. In the simulations, the interesting 
change in coverage fraction occurs at late times, while 
in the experiments, only relatively early times are con- 
sidered 22 1. Also, the fact that the experiments agree 
reasonably well with equilibrium calculations p3| is not 
surprising because the late-time changes in p are rela- 
tively small in the simulations. Therefore, longer experi- 
ments with more precise measurements of p are required 
to directly observe the effects predicted by our simula- 
tions. This discussion naturally leads to the second point: 
how do small changes in coverage fraction produce large 
changes in viscosity? An ad hoc model that is capable of 
explaining the large viscosity rise assumes that the vis- 
cosity is proportional to 1/(A - p), where A is a constant 
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determined by the equilibrium coverage fraction. This 
model is based on the idea that the fluidity (the inverse of 
the viscosity) is proportional to the number of unbound 
sites. Clearly, both more careful direct measurements of 
the coverage fraction versus time and a better theoret- 
ical understanding of the connection between viscosity 
and coverage fraction are needed. 
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FIG. 9. In this figure, the solid points are taken from Fig. 2 
in Ref . j^] . They are the viscosity values (left-hand axis) ver- 
sus time (bottom axis) for pH = 5.5 and a Ca ++ concentration 
in the subphase of 0.65 mM. The solid curve is the simulation 
data for Model A and K — 500. Plotted here is coverage 
fraction p (right-hand axis) versus number of iterations (top 
axis) . 

The final two comments concern possible refinements 
of the adsorption model when applying it to the mono- 
layer system. In this paper, we considered the two cases 
of binding to any open pair of sites (Model A) and bind- 
ing with nearest-neighbor exclusions (Model B) because 
they are simple cases with different geometric arrange- 
ments. The correct detailed description of the Langmuir 
monolayer system is certainly more complicated than ei- 
ther of these. However, as mentioned, the degree of ion- 
ization of the monolayer is pH dependent. To zeroth 
order, Model B is a reasonable description of a mono- 
layer that is only partially ionized for two reasons. First, 
for a partially ionized monolayer, if a particular site is 
available for binding, it is highly unlikely that any of the 
neighbors will be available as well. Second, the steady- 
state values of p found in Model B are in reasonable 
agreement with measurements of the values of p reported 
for monolayers for pH between 5 and 6 p3| . 

The second refinement concerns lateral diffusion of 
Ca ++ ions once they have bound to the monolayer. Inclu- 
sion of diffusion should not substantially alter the quali- 
tative results presented here, but it would effect the quan- 
titative interpretation of the rate constants K + and K~ . 



One can model lateral diffusion of Ca ++ as the unbinding 
of a Ca ++ from one of the monolayer molecules followed 
by a rotation around its remaining bond and subsequent 
binding to another available site. However, this process 
could also be viewed as a complete unbinding and rebind- 
ing at a neighboring site with a renormalized rate con- 
stant. In addition, for p to evolve in time, diffusion would 
need to be coupled with additional binding. This would 
result in rearrangements that are completely analogous 
to transitions from one bad parker to two good parkers. 
Therefore, even with diffusion in the plane of the mono- 
layer, the basic physics remains the same. Jamming will 
still occur, and the slow relaxation of the bad parkers 
due to cooperative behavior will result in the slow time 
scales. 



VI. DISCUSSION 

Equation 2 provides a means of expanding the dynam- 
ics in terms of collective events that occur on slower and 
slower time scales. We were able to directly confirm this 
in the simple situation of the one-dimensional model and 
for the two-particle transitions in Model A. For the one- 
dimensional case, two-particle events were sufficient to 
describe the dynamics of the system, as was found in the 
continuous model. This results in two plateaus in the 
time evolution of the coverage fraction. Single particle 
events, dominated by adsorption, rapidly drive the sys- 
tem to the jamming limit. Processes involving two par- 
ticles are sufficiently slow that p plateaus for some time. 
The length of this plateau is controlled by K, as K ulti- 
mately determines the rates of multiparticle transitions. 
The larger the value of K, the longer the system remains 
at the jamming limit. After enough time, the two par- 
ticle processes have a sufficiently large contribution to 
the dynamics that p increases at a noticeable rate until 
the true steady-state value is reached, and the coverage 
plateaus again. 

In contrast, one can imagine more complicated dynam- 
ics, such as multiple plateaus in the time evolution, occur- 
ring when collective events involving three or more par- 
ticles are important. For example, Fig. 4d illustrates the 
existence of spatial arrangements of unbound sites that 
can not be corrected by two-particle events. In Model 
A, Fig. 8 shows that the two particle events are not ca- 
pable of bringing the system to its steady-state value, as 
they are no longer contributing to the dynamics at late 
enough times. This suggests that the remaining unbound 
sites occur in spatial arrangements that are analogous to 
those in Fig. 4d. Multiple plateaus would arise in the ex- 
treme case where the transition rates for two-particle and 
three-particle events are sufficiently different. This would 
occur as follows. The two-particle transitions would drive 
the system to some value pi in a given time t. If t was 
small enough compared to the three-particle transition 
rate, the system would stay at pi until the three-particle 



7 



events contributed to the dynamics. 

Identifying the existence of multiple plateaus is ex- 
tremely challenging. First, the steady state value of p 
must be sufficiently large that at late times the unbound 
sites are arranged in such a way that two-particle events 
are ineffective. This implies a sufficiently high value of 
K. However, this in turn both decreases significantly the 
rate of collective events and increases the time to reach 
steady state. For Model A and B, we have indirect evi- 
dence of multiple plateaus. In both cases, the coverage 
fraction for K = 10000 appears to be leveling off at a 
value that is lower than the apparent steady-state values 
for K = 500, in the case of Model A, and K = 200, in the 
case of Model B. In principle, p(oo) should approach one 
(or 0.4 for Model B) as K approaches infinity. Therefore, 
the behavior for K = 10000 suggests the beginning of a 
secondary plateau. Unfortunately, as discussed, the time 
required to achieve steady-state increases with K, and 
we do not have sufficient computing power to determine 
if this is a true intermediate plateau for K = 10000 or if 
this is actually the steady-state value. 

It is clear that both analytic and more numeric work 
is needed to fully explore the effects of higher order tran- 
sitions. Identification of the higher order terms in Eq. 2 
is an important step in this process. An exhaustive iden- 
tification of all possible transitions is beyond the scope 
of this paper; however, Fig. 10 identifies a small subset 
of transitions that illustrates why one would expect dif- 
ferences between Model A and B for large enough times 
or large enough values of K. 

Figure 10a shows a set of transitions for Model A, and 
Fig. 10b illustrates the equivalent ones for Model B. In 
both cases, there exists at least two different classes of 
transitions that turn two good parkers (labeled A, B, and 
C in Fig. 10) into one bad parker. For Model A, if A and 
C desorb, then there are two possible sites that result in 
a bad parker, and two possible sites that result in the re- 
establishment of a pair of good parkers. But, if A and B 
desorb, then the situation reverts to the one-dimensional 
case. (In one-dimension, after two good parkers desorb, 
bonding to one out of the three open sites corresponds 
to the creation of a bad parker (see Fig. 4c).) For Model 
B, if A and C desorb, then there are six possible sites 
that result in a bad parker, and six possible sites that 
result in the re-establishment of a pair of good park- 
ers. This results in the same probabilities as in Model 
A. However, in the A to B case, two sites are available 
for bad parkers, and two sites are available for good park- 
ers. Therefore, the chance of two good switching to one 
bad is increased. Because differences in transition rates 
may affect the length of any additional plateaus, detailed 
calculations of these rates are needed for a fuller under- 
standing of the possible dynamics. 
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FIG. 10. This figure illustrates collective events in the 
two-dimensional models, (a) Initially, there are particles at 
the locations labeled A, B, and C. There are two possibilities. 
If A and C desorb, the open circles represent the now avail- 
able sites. If A and B desorb, the open circles at A and the 
gray circles are now the available bonding sites, (b) Again, 
particles are initially at the location A, B, and C. If A and 
C desorb, the open circles represent the now available sites. 
If A and B desorb, the open circles at A and the gray cir- 
cles are now the available bonding sites. In this case, there 
is an additional available site because of the nearest-neighbor 
exclusion. 

In conclusion, we present results of simulations of the 
reversible parking lot model for three different lattices. 
We have directly confirmed the importance of multiparti- 
cle transitions for governing the late time behavior in two 
of the models. The behavior of the third model is consis- 
tent with the other two. We discussed the implications 
of a description of the dynamics in terms of collective 
events. For the right ratios of transition rates, one would 
expect to observe multiple plateaus. There is a sugges- 
tion of intermediate plateaus in our system, but compu- 
tational limits prevented any conclusive evidence. One 
alternative method for finding multiple plateaus would 
be to consider different particle shapes as a means of 
adjusting the relative rates of multiparticle transitions. 
Finally, we presented the possible relevance of the model 
to the binding of Ca ++ to Langmuir monolayers. We 
showed that the jamming and subsequent slow relaxation 
of the binding of Ca ++ ions is a strong candidate for the 
source of the long-time scales observed in the viscosity 
measurements. There are experimental and theoretical 
details that require further exploration, including direct 
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measurements of the Ca ++ coverage fraction, modeling of 
the dependence of viscosity on Ca ++ coverage fraction, 
better modeling of pH effects, and both measurements 
and modeling of lateral diffusion. However, given how 
well the model presented here captures the time scales 
present in the viscosity data, such future studies should 
prove extremely fruitful. 
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